home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / ASTRNOMY / DE118I.ZIP / SQRTL.C < prev    next >
C/C++ Source or Header  |  1992-07-15  |  3KB  |  170 lines
  1. /*                            sqrtl.c
  2.  *
  3.  *    Square root, long double precision
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * long double x, y, sqrtl();
  10.  *
  11.  * y = sqrtl( x );
  12.  *
  13.  *
  14.  *
  15.  * DESCRIPTION:
  16.  *
  17.  * Returns the square root of x.
  18.  *
  19.  * Range reduction involves isolating the power of two of the
  20.  * argument and using a polynomial approximation to obtain
  21.  * a rough value for the square root.  Then Heron's iteration
  22.  * is used three times to converge to an accurate value.
  23.  *
  24.  * Note, some arithmetic coprocessors such as the 8087 and
  25.  * 68881 produce correctly rounded square roots, which this
  26.  * routine will not.
  27.  *
  28.  * ACCURACY:
  29.  *
  30.  *
  31.  *                      Relative error:
  32.  * arithmetic   domain     # trials      peak         rms
  33.  *    IEEE      0,10        30000       8.1e-20     3.1e-20
  34.  *
  35.  *
  36.  * ERROR MESSAGES:
  37.  *
  38.  *   message         condition      value returned
  39.  * sqrt domain        x < 0            0.0
  40.  *
  41.  */
  42.  
  43. /*
  44. Cephes Math Library Release 2.2:  December, 1990
  45. Copyright 1984, 1990 by Stephen L. Moshier
  46. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  47. */
  48.  
  49.  
  50. #include "mconf.h"
  51.  
  52. #define SQRT2 1.4142135623730950488017E0L
  53.  
  54.  
  55.  
  56. long double sqrtl(x)
  57. long double x;
  58. {
  59. int e;
  60. long double z, w;
  61. #ifndef UNK
  62. short *q;
  63. #endif
  64. long double frexpl(), ldexpl();
  65.  
  66. if( x <= 0.0 )
  67.     {
  68.     if( x < 0.0 )
  69.         mtherr( "sqrtl", DOMAIN );
  70.     return( 0.0 );
  71.     }
  72. w = x;
  73. /* separate exponent and significand */
  74. #ifdef UNK
  75. z = frexpl( x, &e );
  76. #endif
  77.  
  78. /* Note, frexp and ldexp are used in order to
  79.  * handle denormal numbers properly.
  80.  */
  81. #ifdef IBMPC
  82. z = frexpl( x, &e );
  83. q = (short *)&x; /* point to the exponent word */
  84. q += 4;
  85. /*
  86. e = ((*q >> 4) & 0x0fff) - 0x3fe;
  87. *q &= 0x000f;
  88. *q |= 0x3fe0;
  89. z = x;
  90. */
  91. #endif
  92. #ifdef MIEEE
  93. z = frexpl( x, &e );
  94. q = (short *)&x;
  95. /*
  96. e = ((*q >> 4) & 0x0fff) - 0x3fe;
  97. *q &= 0x000f;
  98. *q |= 0x3fe0;
  99. z = x;
  100. */
  101. #endif
  102.  
  103. /* approximate square root of number between 0.5 and 1
  104.  * relative error of linear approximation = 7.47e-3
  105.  */
  106. /*
  107. x = 0.4173075996388649989089L + 0.59016206709064458299663L * z;
  108. */
  109.  
  110. /* quadratic approximation, relative error 6.45e-4 */
  111. x = ( -0.20440583154734771959904L  * z
  112.      + 0.89019407351052789754347L) * z
  113.      + 0.31356706742295303132394L;
  114.  
  115. /* adjust for odd powers of 2 */
  116. if( (e & 1) != 0 )
  117.     x *= SQRT2;
  118.  
  119. /* re-insert exponent */
  120. #ifdef UNK
  121. x = ldexpl( x, (e >> 1) );
  122. #endif
  123. #ifdef IBMPC
  124. x = ldexpl( x, (e >> 1) );
  125. /*
  126. *q += ((e >>1) & 0x7ff) << 4;
  127. *q &= 077777;
  128. */
  129. #endif
  130. #ifdef MIEEE
  131. x = ldexpl( x, (e >> 1) );
  132. /*
  133. *q += ((e >>1) & 0x7ff) << 4;
  134. *q &= 077777;
  135. */
  136. #endif
  137.  
  138. /* Newton iterations: */
  139. #ifdef UNK
  140. x += w/x;
  141. x = ldexpl( x, -1 );    /* divide by 2 */
  142. x += w/x;
  143. x = ldexpl( x, -1 );
  144. x += w/x;
  145. x = ldexpl( x, -1 );
  146. #endif
  147.  
  148. /* Note, assume the square root cannot be denormal,
  149.  * so it is safe to use integer exponent operations here.
  150.  */
  151. #ifdef IBMPC
  152. x += w/x;
  153. *q -= 1;
  154. x += w/x;
  155. *q -= 1;
  156. x += w/x;
  157. *q -= 1;
  158. #endif
  159. #ifdef MIEEE
  160. x += w/x;
  161. *q -= 1;
  162. x += w/x;
  163. *q -= 1;
  164. x += w/x;
  165. *q -= 1;
  166. #endif
  167.  
  168. return(x);
  169. }
  170.